
function [lp,glp]=log_posterior_thetas(pars,D,M,X,Z,R,a0,om_a,th0,om_th,b0A,b0D,om_b,sv,dv,f0,om_f,svs,dvs)

[N,T,K]=size(X);
KZ=size(Z,3);
Kth=size(R,2);

alpha=pars(1:N,1);
theta=pars(N+(1:2*Kth),1);
beta=pars(N+2*Kth+(1:2*N),1);
v=pars(3*N+2*Kth+(1:N),1);
f=pars(4*N+2*Kth+(1:N),1);
vs=pars(5*N+2*Kth+(1:N),1);
% xi=pars(6*N+2*Kth+(1:K),1); uniform prior
% gamma=pars(6*N+2*Kth+K+(1:KZ),1); uniform prior
M=reshape(M,N*T,1);
D=reshape(D,N*T,1);
kappa=M.*pars((6*N+2*Kth+K+KZ+1):end,1)./repmat(vs,T,1);

% log_prior
lp=sum(-0.5*((alpha-a0)/om_a).^2)+sum(-0.5*((theta-th0)/om_th).^2)+sum(-0.5*((beta(1:N)-b0A)/om_b).^2)+...
    sum(-0.5*((beta(N+(1:N))-b0D)/om_b).^2)+sum(-(sv+1)*log(v)-dv*v.^(-1))+...
    sum(-0.5*((f-f0)/om_f).^2)+sum(-(svs+1)*log(vs)-dvs*vs.^(-1));

% log_likelihood (use 2nd-order polynomial approx. for tails to ensure finite values)
rl=normpdf(-37.5)/normcdf(-37.5);
ru=normpdf(37.5)/normcdf(37.5,'upper');
lp=lp+sum(log(normcdf(max(-37.5,kappa)).^(M.*D)).*(kappa>=-37.5)+...
    (kappa<-37.5).*(M.*D).*(log(normcdf(-37.5))+rl*(kappa+37.5)-...
    0.5*(-37.5*rl+rl^2)*(kappa+37.5).^2)+...
    log(normcdf(min(37.5,kappa),'upper').^(M.*(1-D))).*(kappa<=37.5)+...
    (kappa>37.5).*(M.*(1-D)).*(log(normcdf(37.5,'upper'))-ru*(kappa-37.5)+...
    0.5*(37.5*ru-ru^2)*(kappa-37.5).^2));

lp=-lp;

if nargout==2
    glp=zeros(length(pars),1);
    % alpha
    glp(1:N,1)=-((alpha-a0)/om_a^2);
    % theta
    glp(N+(1:2*Kth),1)=-((theta-th0)/om_th^2);
    % beta
    glp(N+2*Kth+(1:N),1)=-((beta(1:N)-b0A)/om_b^2);
    glp(2*N+2*Kth+(1:N),1)=-((beta(N+(1:N))-b0D)/om_b^2);
    % v
    glp(3*N+2*Kth+(1:N),1)=-(sv+1)*v.^(-1)+dv*v.^(-2);
    % f
    glp(4*N+2*Kth+(1:N),1)=-((f-f0)/om_f^2);
    % for vs and kappa:
    A=M.*D.*((normpdf(max(-37.5,kappa))./normcdf(max(-37.5,kappa))).*(kappa>=-37.5)+...
        (kappa<-37.5).*(rl-...
        (-37.5*rl+rl^2)*(kappa+37.5)))+...
        M.*(1-D).*((-normpdf(min(37.5,kappa))./normcdf(min(37.5,kappa),'upper')).*(kappa<=37.5)+...
        (kappa>37.5).*(-ru+...
        (37.5*ru-ru^2)*(kappa-37.5)));
    % vs
    glp(5*N+2*Kth+(1:N),1)=-(svs+1)*vs.^(-1)+dvs*vs.^(-2)+...
        sum(reshape(-(kappa./repmat(vs,T,1)).*A,N,T),2);
    % kappa
    glp((6*N+2*Kth+K+KZ+1):end,1)=repmat(vs.^(-1),T,1).*A;
    
    glp=-glp;
end

